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ABSTRACT 

We study the stellar population properties of the IRAC-detected 6 < z < 10 galaxy candidates from the Spitzer 
UltRa Faint SUrvey Program (SURFS UP). Using the Lyman Break selection technique, we find a total of 17 
galaxy candidates at 6 < z < 10 from HST images (including the full-depth images from the Hubble Frontier 
Fields program for MACS1149 and MACS0717) that have detections at S/N > 3 in at least one of the IRAC 
3.6/im and 4.5/im channels. According to the best mass models available for the surveyed galaxy clusters, 
these IRAC-detected galaxy candidates are magnified by factors of ^ 1.2-5.5. Due to the magnification of the 
foreground galaxy clusters, the rest-frame UV absolute magnitudes Mi6oo are between -21.2 and -18.9 mag, 
while their intrinsic stellar masses are between 2 x 10^ Mq and 2.9 x 10 ^Mq. We identify two Lye emitters 
in our sample from the Keck DEIMOS spectra, one at zpya = 6.76 (in RXJ1347) and one at zpya = 6.32 (in 
MACS0454). We find that 4 out of 17 z > 6 galaxy canidates are favored by z ^ 1 solutions when IRAC 
fiuxes are included in photometric redshift fitting. We also show that IRAC [3.6] - [4.5] color, when combined 
with photometric redshift, can be used to identify galaxies likely with strong nebular emission lines or have 
obscured AGN contributions within certain redshift windows. 

Subject headings: galaxies: evolution — galaxies: high-redshift — method: data analysis — gravitational 
lensing 


1. introduction 

Galaxies at 6 < z ^ 10 are one of the frontiers in observa¬ 
tional astronomy because they are a key player in the reioniza¬ 
tion process. It is widely postulated that galaxies provided the 
bulk of ionization photons, but low-level AGN activity (with 
their likely very high escape fractions of ionizing photons) is 
still a possibility (e.g., Giallongo et al. 2015). To improve 
our knowledge about the importance of galaxies on reioniza¬ 
tion, we should measure their ionizing photon production rate 
(through their star formation rate density) and their ionizing 
photon escape fraction (Robertson et al. 2010). We should 
also measure their stellar mass and, under reasonable assump¬ 
tions about their star formation history, infer how many ioniz- 
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ing photons they produced in the past. 

Rest-frame optical stellar emission from galaxies is 
crucial for stellar mass measurement; the rest-frame 4000A 
break shifts to > 3/im in the observed frame and requires 
deep Spitzer observations at the moment. Over a thousand 
z ^ 6 galaxy candidates have been identified in deep HST 
extragalactic blank fields like CANDELS (Koekemoer et al. 
2011; Grogin et al. 2011), HUDF/XDF (Beckwith et al. 2006; 
Koekemoer et al. 2013; Illingworth et al. 2013; Bouwens 
et al. 2015), BoRG (Trend et al. 2011, 2012; Bradley et al. 
2012), and HIPPIES (Yan et al. 2011). Among all the z > 6 
galaxy candidates, more than 100 of them have individual 
Spitzerf[RAC detections (Eyles et al. 2007; Yan et al. 2006; 
Stark et al. 2009; Labbe et al. 2010; Capak et al. 2011; Labbe 
et al. 2013; Roberts-Borsani et al. 2015; Labbe et al. 2015); 
their IRAC fiuxes enable more robust constraints on their stel¬ 
lar masses. The inferred stellar masses of these z ^ 6 galaxy 
candidates range from ^ 10^ to ^ 10^^ Mq, surprisingly large 
for a universe younger than 1 Gyr old. It is likely that the ma¬ 
jority of IRAC-detected galaxies are at the high-mass 
end of the stellar mass function, although some of these galax¬ 
ies likely have their IRAC fiuxes boosted by strong nebular 
emission lines like [OIII], Hg, and H/3 (e.g., Finkelstein et al. 
2013; Smit et al. 2014; De Barros et al. 2014). 

Using the strong gravitational lensing power of rich galaxy 
clusters is a novel avenue to explore high-redshift galaxies 
(Soucail 1990). Galaxy candidates at z ^ 6 that are magni¬ 
fied by foreground clusters were starting to be identified from 
HST images more than a decade ago (e.g., Ellis et al. 2001; 
Hu et al. 2002; Kneib et al. 2004); these observations pro¬ 
vide an alternative way to probe the faint-end of the luminos¬ 
ity function with shorter exposure time than in blank fields. 
Recently, the Cluster Lensing And Supernova survey with 
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Hubble (CLASH; Postman et al. 2012) program and HST- 
GO-11591 (PI: Kneib) program observed 34 galaxy clusters. 
The ongoing Hubble Frontier Fields (HFF; PI: Lotz^^) pro¬ 
gram, upon its complete execution, will obtain deep HST 
ACSAVFC3-IR images in six galaxy cluster fields (four of 
them are in the CLASH sample). Bradley et al. (2014) re¬ 
cently reported 262 6 < z < 8 galaxy candidates across 18 
clusters in the CLASH sample based on photometric red- 
shift selection, and they demonstrated the power of using 
strong gravitational lensing to identify high-z galaxies, espe¬ 
cially at the bright end of the luminosity function. The even 
deeper HFF data, despite being >0.7 mag shallower than the 
HUDF data, can probe galaxies intrinsically fainter than can 
be probed in the HUDF due to the power of gravitational lens¬ 
ing. Other coordinated campaigns are also underway to com¬ 
plement the deep HST images in those targeted galaxy clus¬ 
ter fields (e.g., the Grism Lens-Amplified Survey from Space 
program; Schmidt et al. 2014; Treu et al. 2015). 

Here we use the deep SpitzerHRAC images obtained from 
the Spitzer UltRa Faint SUrvey Program (SURFS UP; Bradac 
et al. 2014; hereafter Paper I) to probe the rest-frame optical 
emission from z > 6 galaxy candidates. SURFS UP surveys 
10 strong-lensing cluster fields (our sample partially overlaps 
with both CLASH and HFF) with Spitzer IRAC images in 
the 3.6pm and 4.5pm channels, with exposure times of > 28 
hours per channel per cluster. Paper I summarizes the sci¬ 
ence motivations and observational strategies of SURFS UP, 
and Ryan et al. (2014) presented the z ~ 7 galaxy candidates 
in the Bullet Cluster, one of which is detected in both IRAC 
channels. In this work, we explore the 6 < z ^ 10 galaxy can¬ 
didates with IRAC detections in 8 additional clusters in our 
sample^^ and present their physical properties inferred from 
their broadband fiuxes. We also make all the IRAC imaging 
data available for the community on our webpage 

The structure of this paper is as follows: Section 2 describes 
our HST and Spitzer imaging data and photometry; Section 
3 describes our 6 < z ^ 10 galaxy sample selection proce¬ 
dure; Section 4 describes the identification of Tya emission 
from spectroscopy for two galaxy candidates with IRAC de¬ 
tection at z = 6.76 (RXJ1347-1216) and z = 6.32 (MACS0454- 
1251); Section 5 presents our spectral energy distribution 
(SED) modeling procedure and results, and Section 6 explores 
the idea of using IRAC [3.6] - [4.5] color to identify galaxies 
with strong nebular emission lines. Finally, Section 7 sum¬ 
marizes our findings. Throughout the paper, we assume a 
ACDM concordance cosmology with = 0.3, Ua = 0.7, and 
the Hubble constant Hq = 70kms“^ Mpc“^ Coordinates are 
given for the epoch J2000.0, and all magnitudes are in the AB 
system. 

2. IMAGING DATA AND PHOTOMETRY 
2.1. HST Data and Photometry 

We list the eight galaxy clusters analyzed in this work 
in Table 1. Among the eight clusters, six (MACS0717, 
MACS0744, MACS1149, RXJ1347, MACS1423, and 
MACS2129) are in the Cluster Lensing And Supernova Sur¬ 
vey (CLASH; Postman et al. 2012) sample; therefore, each of 

http://www.stsci.edu/hst/campaigns/ 

frontier-fields 

The HST WFC3/IR imaging data for the tenth cluster, MACS2214, will 
be obtained in late 2015 {HST-GO 13666; PI: Bradac). 

http : / /vivivf. physics . ucdavis . edu/aLi jmarusa/ 
SurfsUp.html 


them has HST imaging data in at least twelve ACSAVFC and 
WFC3/IR filters^^ from the ACS and WFC3/IR cameras. The 
typical 5(7 depths for the CLASH clusters reported by Post¬ 
man et al. (2012) are between 27.0 and 27.5 mag (within 0^'4 
diameter apertures) for each filter, and the large number of 
filters provides unique constraints for high-z galaxy searches 
among the HST deep fields. 

In addition to the CLASH data, full-depth images for 
MACS0717 and MACS1149 from the Hubble Frontier Fields 
(HFF) program have also been released in June 2015. In these 
two clusters, HST spent a total ^140 orbits that are roughly 
split between ACS and WFC3/IR filters, and these images 
achieve ^ 28.7-29 mag^^ in the optical (ACS) and NIR 
(WFC3). We use the deepest HFF images for MACS 1149 
and MACS0717 for photometry in the filters where such im¬ 
ages are available^^. The six CLASH clusters in our sam¬ 
ple are also observed by the Grism Lens-Amplified Survey 
from Space (GLASS; PI: Treu; Schmidt et al. 2014; Treu 
et al. 2015) program (MACS0717, MACS1149, MACS1423, 
RXJ1347, MACS0744, and MACS2129) and have deep grism 
spectra available. 

For the remaining two clusters being analyzed in this work, 
we have data for RCS2327 as part of the SURFS UP HST 
observations (F^5'r-G013177 PI Bradac; Hoag et al. 2015) 
and previous archival data (F^5T-GO 10846 PI Gladders; see 
also Sharon et al. 2015). For MACS0454, we use the archival 
observations from i75r-G011591 (PI: Kneib), GO-9836 (PI: 
Ellis), and GO-9722 (PI: Ebeling). We list the limiting mag¬ 
nitudes of point sources within a 0."4 aperture for MACS0454 
and RCS2327 in Table 2. 

We will use the template-fitting software T-PHOT (Merlin 
et al. 2015) — the successor to TEIT (Laidler et al. 2007) — 
to measure the colors between HST and Spitzer IRAC images 
(see Section 2.2). To prepare the HST images for T-PHOT, 
we use the public 0."03/pixel CLASH images and block-sum 
the images to make 0."06/pixel images. We also edit the astro- 
metric image header values (CRVALs and CRPIXs) to con¬ 
form to T-phot’s astrometric requirements^^ and make sure 
that HST and Spitzer images are aligned to well within 0^'l 
(Paper I). 

We also match the point-spread functions (PSPs) among all 
HST filters to get consistent colors. To do so, we identify 
isolated point sources in each cluster field, and we use the 
p s f mat ch task in I RAF to match all HST images to have the 
same PSP as the reddest band, P160W. In practice, because 
of the small field of view of each cluster and the crowded 
environment, we can only select ^ 5 isolated point sources 
in each cluster for PSP-matching. However, we measure the 
curves of growth of each point source after PSP-matching and 
find that in most filters, the curves match to within ^20% of 
thatofP160W. 

After pre-processing, we extract photometry on HST im¬ 
ages using SExtractor (Bertin & Arnouts 1996; version 2.8.6). 
We use the combined IR images as the detection image, 
and the SExtractor detection/deblending settings similar to 

For the CLASH clusters, the ACS filters include F435W, F475W, 
F606W, F625W, F775W, F814W, and F850LP; the WFC3/IR filters include 
F105W, FI low, F125W, F140W, and F160W. 

1^ The magnitude limits are from the Frontier Fields website. 

1^ The HFF filter sets include F435W, F606W, F814W, F105W, F125W, 
F140W, and F160W. 

18 T-PHOT requires that the pixel boundaries of the high- and low- 
resolution images be perfectly aligned, meaning that both images should have 
the same CRVALs and both have half-integer CRPIXs. 
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TABLE 1 

SLFRFS UP Galaxy Cluster Sample 



Cluster Name 

Short Name^ 

R.A. 

(deg.) 

Decl. 

(deg.) 

-^cluster^ 

A^bg" 

A^lbgjrac^ 

1 

MACS J0454.1-0300 

MACS0454 

73.545417 

-3.018611 

0.54 

10 

2 

2 

MACSJ0717.5+3745^’f 

MACS0717 

109.390833 

37.755556 

0.55 

10 

0 

3 

MACSJ0744.8+3927^ 

MACS0744 

116.215833 

39.459167 

0.70 

4 

1 

4 

MACSJ1149.5+2223^’f 

MACS 1149 

177.392917 

22.395000 

0.54 

11 

3 

5 

RXJ1347-1145^ 

RXJ1347 

206.883333 

-11.761667 

0.59 

9 

3 

6 

MACSJ1423.8+2404^ 

MACS 1423 

215.951250 

24.079722 

0.54 

9 

6 

7 

MACSJ2129.4-074U 

MACS2129 

322.359208 

-7.690611 

0.59 

0 

0 

8 

RCS2-2327.4-0204 

RCS2327 

351.867500 

-2.073611 

0.70 

6 

1 

9 

1E0657-56 

Bullet Cluster 

104.614167 

-55.946389 

0.30 

10 

1 

10 

MACS2214.9-1359S 

MACS2214 

333.739208 

-14.003000 

0.50 

N/A 

N/A 

Total 
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^ We will refer to each cluster by its short name. 

^ Cluster redshift 

Number of 6 < z < 10 LBG candidates selected by their HST colors. 

^ Number of 6 < z < 10 LBG candidates with > 3a detections in at least one IRAC channel. 
^ A CLASH cluster 
^ A Hubble Frontier Fields cluster 

§ The HST WFC3/IR data for MACS2214 will be collected in late 2015. 


TABLE 2 

HST 5a Limiting magnitudes (point source, Of'4 aperture) eor RCS2327 and MACS0454 


Cluster E435W F555W F775W F814W F850LP F098M FllOW F125W F160W 


MACS0454 ••• 27.7 26.9 27.9 26.5 ••• 28.1 ••• 27.4 

RCS2327 27.6 • • • • • • 27.6 • • • 27.3 • • • 27.6 27.5 


(but slightly more conservative than) the values adopted by 
CLASH for their high-z galaxy search (Postman et al. 2012). 
Because our focus in this work is on identifying IRAC- 
detected high-z sources, the slightly more conservative set¬ 
tings do not reject many potential IRAC-detected candidates 
but eliminates most spurious detections. We also follow the 
procedure outlined in Trenti et al. (2011) to rescale the flux 
errors reported by SExtractor. At the end of the process, we 
use the resulting photometric catalogs and segmentation maps 
for both IRAC photometry (Section 2.2) and color selection 
(Section 3). 

2.2. IRAC Data and Photometry 

The IRAC imaging data set for SURFS UP was presented 
in Paper I; the total exposure time for each IRAC channel is 
about 30 hours (see Ryan et al. 2014 and Paper I for details.) 
The coadded IRAC mosaics are deeper within the main cluster 
helds covered by WFC3/IR, and the typical 3cr limiting mag¬ 
nitude within 3”-radius apertures is 26.6 mag in the 3.6 jim 
channel (hereafter chi) and 26.2 mag in the 4.5 jim channel 
(hereafter ch2) where source blending is not severe. 

We use T-PHOT to measure consistent colors between HST 
and Spitzer IRAC images with a template-htting approach 
(see also Laidler et al. 2007 for the template-htting concept 
employed by T-PHOT.) The template-htting approach has 
been demonstrated to work well for blank-held surveys such 
as CANDELS (Guo et al. 2013), but it does require images 
with zero mean background. Most galaxy cluster helds have 
considerable spatial variations in local sky background, so 
subtracting a constant background does not generally work. 
Therefore, instead of htting ah sources in the held at the same 
time — which is the strategy for blank-held surveys — we 
subtract the local background and perform the ht for each 
high-z candidate separately to get the cleanest residual pos¬ 


sible. 

As described in Merlin et al. (2015), T-PHOT is designed 
to measure the huxes in the low-resolution image (in our case 
the IRAC images) for ah the sources detected in the high- 
resolution image (in our case the F160W images). T-PHOT 
does so by constructing a template for each source; it con¬ 
volves the cutout of each source in the F160W image by a 
PSF-transformation kernel that matches the F160W resolution 
to the IRAC resolution. Once the templates are available (and 
with their huxes normalized to 1), T-PHOT solves the set 
of linear equations and hnds the combination of coefficients 
for each template that most closely reproduce the pixel values 
in IRAC images; each coefficient is therefore the hux of the 
source in IRAC. T-PHOT also calculates the full covariance 
matrix and uses the diagonal terms of the covariance matrix 
to calculate hux errors. For each source, T-PHOT reports 
a “covariance index”, dehned as the ratio between the max¬ 
imum covariance of the source with its neighbors (max(cr/y)) 
over its hux variance {(Tu), which serves as an indicator of 
how strongly correlated the source’s hux is with its closest 
(or brightest) neighbor. Generally, a high covariance index 
(> 1) is associated with more severe blending and large hux 
errors, at least from simulations (Laidler et al. 2007; Merlin 
et al. 2015). Therefore, sources with high covariance indices 
should be treated with caution. 

Obviously the PSF-transformation kernel that matches the 
F160W PSF to the IRAC PSF is a crucial element in this 
process. We generate IRAC PSFs by stacking point sources 
observed in the exposures from both the primary cluster 
held and the hanking held. We identify point sources using 
Sextractor with DEBLEND_MINCONT=10“'', MINAREA=9, 
DETECT_THRESH=ANALYSIS_THRESH=2, and a Gaus- 
sian convolution kernel with (7 = 3 pixels (dehned over a 5 
by 5-pixel grid). We require that ah point sources have an 
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Fig. 1.— Brightness and half-light radii for all sources in IRAC chi of the 
SURFS UP clusters. The half-light radii (FLUX_RADIUS) are in 0'.'6 IRAC 
pixels. The box illustrates the crudely defined stellar loci where all point 
sources are expected to fall. When selecting putative point sources for each 
field, we first make sure the objects are near this locus for a given cluster, 
then place additional constraints on proximity to neighbors, axis ratio, and 
re-centering/alignment accuracy for the final sample per cluster. There are 
typically > 40 stars per cluster for PSF generation. 

axis ratio oihja> 0.9, lie on the stellar locus within the box 
shown in Figure 1 , and are sufficiently separated from neigh¬ 
boring objects to have reliable centroids (FLAGS< 3). We 
recompute the PSF centroids by fitting a Gaussian profile to 
the inner profile (r < 4 pixels) using the Sextractor barycen- 
ters as initial guesses, and align the point sources using sine 
interpolation. To mask neighboring objects, we grow the seg¬ 
mentation maps from Sextractor by 2 pixels. At each phase 
we subtract the local sky (assuming there are no local gradi¬ 
ents) and normalize the flux of the point source to unity. We 
sigma-clip average the masked, registered, normalized point 
sources and do one further background correction only to en¬ 
sure the convolutions with T-PHOT are flux conserving. As 
discussed in Paper I, our stacked PSFs are consistent with the 
IRAC handbook. Each of our clusters contains at least 40 
point sources per bandpass in our PSF-making process. 

In practice, T-PHOT still breaks down in very crowded 
regions (e.g., near the cluster center or near bright cluster 
galaxies); in this case, we are limited mostly by our knowl¬ 
edge of the IRAC PSFs and our ability to subtract sky back¬ 
ground underneath the sources. We also measure the “re¬ 
duced for each source in IRAC within a 2f'4 by 2."4 
box by calculating the average difference per pixel between 
the model pixel values and the observed pixel values: xj, = 

Ef‘’“(Amodei-/;',obs)^/(Abs where /;•,model and /;-,obs 

are the model (best-fit) and observed fiux in pixel /, respec¬ 
tively. Later in this work, we only report the IRAC fiuxes 
of the high-z galaxies with reliable IRAC fiux measurements, 
i-e., < 3. 

For the sources with nominal 5/A above 3 but with poor 
T-PHOT residuals, we do not trust the T-PHOT-measured 
fiuxes and estimate the local 3 a fiux limits via artificial source 
simulations. We insert artificial point sources into the F160W 
image (but not into the IRAC image) within 5" of the high-z 
candidates and run T-PHOT to measure the local sky level. 
We repeat this process at least 100 times near each high-z 
candidate and use the resulting IRAC fiux distribution to de¬ 
termine the 3cr fiux limits. In our analysis in Section 5, we use 
the 3(7 fiux limit for only one IRAC filter for one source (chi 
for MACS 1423-1384); for all the other sources, their IRAC 
fiux measurements in both IRAC channels pass the x?. test. 


We also run a separate set of simulations that indepen¬ 
dently estimate the magnitude errors in case T-PHOT under¬ 
estimates the magnitude errors in crowded regions even when 
they pass the xj, test. In this set of simulations, we insert fake 
point sources around each high-z target (within 5") in IRAC 
images with the same magnitude as the T-PHOT-measured 
value, and measure the fiux of the fake sources again with 
T-PHOT. We then calculate the median difference between 
the input and output magnitudes of the fake sources as an in¬ 
dependent magnitude error estimate. We find that for most 
sources, the T-PHOT-reported magnitude error is within 0.1 
mag from the simulated magnitude error, but sometimes the 
simulated magnitude error is much larger than the T-PHOT- 
reported value. In these cases, T-PHOT might underestimate 
the true magnitude errors, so we adopt the simulated mag¬ 
nitude errors in our SED modeling. We note that by adding 
fake sources in IRAC images, we increase the fiux error due to 
crowding, so the simulated magnitude errors could be higher 
than the true magnitude errors. 

3. SAMPLE SELECTION 

We select galaxy candidates at 6 based on their rest- 
frame UV colors using the Lyman-break selection method 
(Steidel & Hamilton 1993; Giavalisco 2002). For the CLASH 
clusters, we use the published color criteria presented below 
for selecting z ^ 6, 7, 8, and 9 Lyman-break galaxies (LBGs); 
for RCS2327 and MACS0454, we design our own color cuts. 
All galaxy colors are calculated using their isophotal mag¬ 
nitudes (MAG_ISO) from SExtractor. After the initial color 
selection, we inspect the galaxy candidates to remove image 
artifacts and objects with problematic photometry. We then 
measure each LBG candidate’s fiuxes in IRAC, and we only 
present the candidates with S/N > 3 in at least one channel. 
Becuase of the differences in the available filters in each clus¬ 
ter, we explain our color-selection process in more detail be¬ 
low and list the sample size in Table 1 . The full sample of the 
color-selected z ^ 6 galaxy candidates with IRAC detections 
is presented in Table 3. 

3.1. CLASH & Hubble Frontier Fields Clusters: MACS0717, 
MACS 1149, MACS0744, MACS1423, MACS2129, and 
RXJ1347 

For the six clusters that in the CLASH sample, we use the 
criteria below. To be selected as a z ^ 6 LBG candidate, a 
source has to satisfy all of the following criteria from Gonza¬ 
lez etal. (2011): 

^ F775W-F850LP> 1.3 
F850LP-F125W<0.8 
S/N>5 inF850LPandF125W 
^ S/N < 2 in filters bluer than F606W 

where we calculate all signal-to-noise ratios (S/N) within 
the isophotal (ISO) aperture^^. If the S/N in F775W is be¬ 
low one, we use the la fiux limit in F775W to calculate the 
F775W-F850LP color. To select LBG candidates at z ^ 7, 
we use the color criteria from Bouwens et al. (201 1): 

The S/N limits in blue filters roughly correspond to magnitude limits 
of > 28 mag, based on the typical limiting magnitudes presented by Postman 
etal. (2012). 
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TABLE 3 

IRAC-Detected 6 < z < 10 Galaxy Candidates 


Object ID 

R.A. 

Deck 

F160W' 

[3.6]2 

«3.6^ 

[4.5]* 

^ 4 . 5 ^ [3.6]-[4.5] Spectroscopy^ 


(deg.) 

(deg.) 

(mag) 

(mag) 


(mag) 

(mag) 


F125W-dropouts (z ~ 9) 


MACS1149-JD^ 

177.389943 

22.412719 

25.6±0.1 

25.8 ±0.4 

0.27 

25.0 ±0.2 

0.29 

0.8±0.4 

M,G 

F105W-dropouts (z ~ 8) 

MACS1423-1384 

215.942115 

24.079401 

25.7 ±0.2 

> 23.6*' 

0.95 

24.1 ±0.4^ 

0.95 

>-0.5 

G 

RXJ1347-1080h 

206.891236 

-11.752594 

26.3 ±0.2 

25.4 ±0.2 

0.12 

25.3 ±0.2 

0.11 

0.2 ±0.2 

D,G 

F850LP-dropouts (z ~ 7) 

MACS0744-2088h 

116.250405 

39.453011 

25.5 ±0.2 

25.2 ±0.4^ 

0.50 

25.1 ±0.2 

0.50 

0.1 ±0.3 

G 

MACS 1423-587 

215.940493 

24.090848 

25.3±0.1 

24.3 ±0.3^ 

0.86 

26.2 ±0.5 

0.80 

-L8±0.6 

G 

MACS 1423-774 

215.935607 

24.086475 

25.9 ±0.2 

25.1 ±0.2 

0.70 

25.5 ±0.3 

0.69 

-0.4 ±0.3 

D,G 

MACS 1423-2248 

215.932958 

24.070875 

25.6±0.1 

25.0±0.1 

0.42 

25.3 ±0.2 

0.45 

-0.2 ±0.2 

D,G 

MACS1423-1494 

215.935871 

24.078414 

26.3 ±0.2 

26.1 ±0.4 

0.92 

25.2 ±0.2 

0.89 

0.9 ±0.4 

D,G 

MACS1423-2097d 

215.945534 

24.072433 

25.8 ±0.2 

24.6 ±0.3^ 

0.68 

24.6 ±0.1 

0.68 

0.0 ±0.2 

D,G 

RXJ1347-1216d’^’h 

206.900848 

-11.754199 

26.1 ±0.2 

24.3 ±0.1 

0.16 

25.6 ±0.2 

0.11 

-L3±0.2 

D,G 

RXJ1347-1800 

206.881657 

-11.761483 

25.4 ±0.2 

24.3 ±0.3 

0.63 

25.6 ±0.7 

0.55 

-L3±0.8 

G 

Bullet-3g 

104.667375 

-55.968067 

25.0 ±0.2 

23.8 ±0.3 

n/a 

23.8 ±0.3 

n/a 

0.0 ±0.4 

FORS2 

F814W-dropouts (z ~ 6-7) 

RCS2327-1282 

351.880595 

-2.076292 

24.8 ±0.1 

24.4 ±0.1 

0.08 

24.1 ±0.1 

0.06 

0.3±0.1 

D,M 

MACS0454-125F 

73.535653 

-3.004116 

24.1 ±0.1 

23.2 ±0.2 

0.60 

23.4 ±0.2 

0.60 

-0.2 ±0.2 

D 

MACS0454-1817 

73.551806 

-3.001018 

26.4 ±0.2 

24.1 ±0.3^ 

0.31 

24.5 ±0.2^ 

0.31 

-0.4 ±0.1 

D 

F775W-dropouts (z - 6) 

MACS 1149-274*' 

177.412009 

22.415783 

24.8 ±0.04 

24.1 ±0.1 

0.44 

24.0±0.1 

0.36 

0.0±0.1 

G 

MACS1149-1204*' 

177.378959 

22.402429 

25.0±0.1 

24.3 ±0.1 

0.64 

24.4 ±0.1 

0.67 

-0.1 ±0.1 

G 


^ Lensed total magnitude (MAG_AUTO) in F160W; the magnification factors (fi) are listed in Table 4. 

^ Isophotal lensed magnitude in IRAC channel 1 based on the isophotal aperture defined in F160W. 

^ /? 3.6 and /? 4.5 are the covariance indices for the [3.6] and [4.5] measurements, respectively. The covariance index of a source i is defined as the ratio between the 
maximum covariance among the neighbors (aij) over the flux variance of itself (an) in the covariance matrix. 

^ Isophotal lensed magnitude in IRAC channel 2 based on the isophotal aperture defined in F160W. 

^ Instruments we used for spectroscopy: D=DEIMOS, M=MOSFIRE, G=HST grism from the GLASS program (Schmidt et al. 2014; Treu et al. 2015, in preparation) 

^ First reported by Zheng et al. (2012). 

^ The IRAC residual in the 3.6^im channel shows that T-PHOT breaks down due to severe blending, so we report the simulated 3cr magnitude limit. 

^ T-PHOT likely underestimates the magnitude errors for these sources due to crowding, so we use the simulated magnitude errors. 

^ Also reported by Smit et al. (2014). 

^ Has Lya detection at z = 6.76. 

^ Has tentative Lya detection at z = 6.32. 

® Reported by Ryan et al. (2014). 

^ Also reported by Bradley et al. (2014). 

criteria from Zheng et al. (2014): 

j^F125W-F160W>0.8 

i 5/A^>5inF160W (4) 

[ < 2 in filters bluer than F850LR 

In total, we find a total of 43 z > 6 LEG candidates from 
the six CLASH/HFF clusters, among them 13 have > 3cr de¬ 
tections in at least one IRAC channel. 

3.2. MACS0454 

For the two galaxy clusters that are not in the CLASH sam¬ 
ple (MACS0454 and RCS2327), we design our own selec¬ 
tion criteria for z ^ 6 galaxy candidates. For MACS0454, we 
have HST imaging data in F555W, F775W, F814W, F850LP, 
FI low, and F160W, although the images in F115W and 
F850LP are shallower than the other filters and we use both 
filters only for S/N rejection of low-z interlopers. We use 
the following criteria to select 6 < z ^ 7.5 galaxy candidates 


^F850LP-F105W>0.7 

F105W-F125W<0.45 

F850LP-F105W > 1.4 x (F105W-F125W) + 0.42 
' S/N>5 inF105WandF125W 
S/N <5 inF814W 
^ S/N <2 in filters bluer than F775W 

To select the LBGs at z ^ 8, we use the criteria from 
Bouwens et al. (2011): 

{ F105W-F125W>0.45 
F125W-F160W<0.5 
5/A^>5inF125WandF160W 
S/N <2 in filters bluer than F850LR 

Finally, to select the LEG candidates at z ^ 9, we use the 

















6 



Fig. 2.— Color-color diagram of the F814W-dropout selection from 
MACS0454. The shaded region (light blue) shows where the expected 
F814W-dropout colors should be. We also plot (in solid curves) the theo¬ 
retical color tracks of a 100 Myr old stellar population with constant star for¬ 
mation taken from Bruzual & Chariot (2003) with different amounts of dust 
attenuation. Sources within the shaded region that also pass the S/N cuts are 
shown in unfilled squares; two of them, MACS0454-1251 and MACS0454- 
1817 (shown in large filled symbols), are detected in IRAC. The color tracks 
of z < 3 galaxies, calculated from the local galaxy templates of Coleman 
et al. (1980), are shown in dashed curves. We also show the expected colors 
of stars from Pickles (1998). 



Fig. 3.— Color-color diagram of the F814W-dropout selection from 
RCS2327. The plot style and model assumptions are the same as Figure 
2. A total of six sources satisfy the F814W-dropout color criteria, and one of 
them (RCS2327-1282; shown as a red circle) has IRAC detections. 


3.3. RCS2327 

For RCS2327, we have deep HST images in F435W, 
F814W, F098M, F125W, and F160W, so we use the following 
criteria for 6 < z < 7.5 galaxy candidates (F814W-dropouts): 

'F814W-F098M>2.2 
F098M-F160W< 1.6 

< F814W-F098M>2.0x(F098M-F160W) + 1.0 (6) 

S/N > 5 in F098M and F160W 
J/N <3 in F435W. 

When the S/N in F814W is less than unity, we use its la 
flux limit to calculate colors. We demonstrate the F814W- 
dropout selection from RCS2327 in Figure 3, and six sources 
pass the above color and signal-to-noise cuts. One of the 
six sources, RCS2327-1218, is detected in both IRAC chan¬ 
nels. In addition to the above criteria, we also use the criteria 
similar to the BoRG survey (Trent! et al. 2011) to search for 
z ^ 7.5 galaxy candidates: 

F098M-F125W> 1.75 
S/N > 5 in F125W and F160W (7) 

S/N < 3 in F814W and F435W. 

and when the S/N in F098M is less than unity, we use its 
Icr flux limit to calculate colors. The F098M-dropout search 
yields no galaxy candidate, so we And a total of one galaxy 
candidate with IRAC detections at z > 6 from RCS2327 (see 
also Hoag et al. 2015 for more details on the dropout search 
in RCS2327). 

To summarize, we And a total of 69 LBG candidates at z ^ 6 
from 9 clusters in SURFS UP; 17 of them have IRAC detec¬ 
tions in at least one channel. Figure 4 shows the cutouts of 
the 16 IRAC-detected LBG candidate in HST and Spitzer im¬ 
ages (one candiate was reported by Ryan et al. 2014). We also 
report the IRAC photometry for the entire sample in Table 3. 
We use the simulated IRAC magnitude errors for MACS 1423- 
1384, MACS1423-587, MACS0744-2088, MACS 1423-2097, 
and MACS0454-1817, because we And that T-PHOT likely 
underestimates their IRAC magnitude errors from the simu¬ 
lations. We keep MACS 1423-1384 in our sample because it 
has a nominal 5.9cr detection in ch2 from T-PHOT, although 
the simulated magnitude error suggests that the additional flux 
error due to crowding reduces it to a 2.2cr detection in ch2. 


(F814W-dropouts): 

'F814W-F110W> 1.0 
F110W-F160W<0.3 
^ 5/V > 5 in FI low and F160W 
^ S/N <2 in F555W. 

We show the HST color-color diagram of the F814W-dropout 
selection in Figure 2. A total of ten sources satisfy the color 
and S/N cuts listed above; among them, MACS0454-1251 and 
MAC0454-1817 are detected in at least one IRAC channel. 

Because we have only one Alter redward of FI lOW, we re¬ 
frain from searching for FI lOW-dropouts in MACS0454 as it 
would yield objects detected in only one HST Alter. 
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Fig. 4.— Cutouts of the first eight IRAC-detected, z > 6 LEG candidates in SURFS UP (excluding the one candidate in the Bullet Cluster reported by Ryan 
et al. 2014). Each row shows the cutout in two HST ACS filters (F435W and F814W in all clusters except for MACS0454, where we show F435W and F555W), 
two HST WFC3/IR filters (F125W and F160W), and two IRAC channels. We also show the neighbor-subtracted cutouts around each LEG candidate in both 
IRAC channels (designated by CH1_RESID and CH2_RESID, respectively). The LEG candidate ID is to the left of each row. Each panel is centered on the 
LEG candidate (marked by the red lines), and each panel spans 10" by 10" on the sky. 
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Fig. 5.— Same as Figure 4, but for the remaining seven LEG candidates. 
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4. SPECTROSCOPIC OBSERVATIONS 

We report the detection of two likely hya emitters among 
our sample with DEIMOS (Faber et al. 2003) on the Keck II 
telescope. The DEIMOS observation is part of a larger cam¬ 
paign to systematically target lensed high-z galaxies behind 
strong-lensing galaxy clusters with DEIMOS and MOSFIRE 
(McLean et al. 2010; McLean et al. 2012) on Keck.We ob¬ 
served six cluster fields between 2013 April and 2014 Au¬ 
gust and targeted 9 out of 17 high-z galaxy candidates in 
Table 3 with DEIMOS. Which galaxy candidates were ob¬ 
served with DEIMOS and MOSFIRE are indicated in Table 3. 
The DEIMOS data were reduced using the standard DEEP2 
spec2d pipeline slightly modified to reduce the data observed 
also with 600 1/mm and 830 1/mm gratings (Lemaux et al. 
2009; Newman et al. 2013). We focus here on the two galax¬ 
ies (RXJ1347-1216 and MACS0454-1251) that have line de¬ 
tections; we will present the full spectroscopic survey (with 
both DEIMOS and MOSFIRE) and the line fiux limits for the 
non-detections in a future work. In addition to the Keck ob¬ 
servations, 13 of the 17 galaxy candidates in Table 3 are also 
observed by the Grism Lens-Amplified Survey from Space 
(GLASS; HST GO-13459; PI: Treu) program; the spectro¬ 
scopic constraints on Lya emission from the HST grism data 
will be presented by Schmidt et al., 2015 (in preparation). 

Below we discuss the two galaxy candidates with robust 
line detections and the likelihood that they are Lya emitters 
at z = 6.76 and z = 6.32. 

4.1. RXJ1347-1216 

We selected this object as a z ^ 7 LBG candidate for 
spectroscopic follow-up on April 3 2013 and May 26 2014. 
This source was also selected by both Smit et al. (2014) 
and Bradley et al. (2014) as a Zphot ^ 6.7 LBG with high 
[OIII]+H/3 equivalent widths (> 1300A rest frame). We used 
the 8301/mm grating in the 2013 run and the 1200 1/mm grat¬ 
ing in the 2014 run, and the total integration times are roughly 
6000 and 7200 seconds, respectively. We had good (but not 
photometric) conditions with ^1" seeing in the 2013 run, but 
the conditions were highly variable in the 2014 run. There¬ 
fore, we only present the line fiux measurements from the 
2013 run, though the line was detected significantly in both 
runs. 

In Figure 6 we show the two-dimensional spectra of 
RXJ1347-1216 from both observation runs (in the top and 
middle panels) and the combined one-dimensional spectrum 
(in the bottom panel). We detect an extended emission feature 
with FWHMobs ^ 16.5 A in the 2013 run, and although the 
blue side of the feature is severely contaminated by sky line 
residual, its asymetric profile with a tail to the red side of the 
spectrum strongly suggests that it is Lya. Using the centroid 
of the sky line residual at 9439A as the peak of of line pro¬ 
file, we determine its Lya redshift as ZLya = 6.76 ± 0.003 (the 
uncertainty corresponds to the width of the sky line residual). 
The Lya redshift is in excellent agreement with its photomet¬ 
ric redshift Zphot = 6.8 ± 0.1, lending additional support to the 
identification of the Lya feature. 

The emission feature is also independently detected at 4cr 
in the GLASS grism data at ^ 9440A. With the grism spectra 
in both G102 and G141, one can test the possibility of this 
feature being an [Oil] doublet (A3727, A3729) at z ^ 1.5 by 
looking for the [OIII] A5007 line at ^ 1.3/im. For a typical 
star forming galaxy, the line ratio [OIII]/[OII] should be at 
least ^ 0.3 (Jones et al. 2015), and the [OIII]/[OII] ratio is 



Aobs [A] 



9400 9420 9440 9460 9480 9500 


Aobs [A] 

^obs [^] 


9400 9420 9440 9460 9480 9500 



^rest [^] 

Fig. 6.— Reduced two-dimensional (top and middle panels) and one¬ 
dimensional (bottom panel) spectrum of RXJl347-1216. The top panel 
shows the data taken with the 830 1/mm grating on April 3, 2013, and the 
middle panel shows the data taken with the 1200 1/mm grating on May 26, 
2014. The one-dimensional spectrum was extracted from the 8301/mm spec¬ 
trum. We also plot the RMS spectra in dashed lines and mark the emission 
line redshift if the line is to be Lya. The flux density values given on the 
ordinate are calculated in the rest-frame assuming the line to be Lya. 

even higher for low-metallicity galaxies (e.g., Maiolino et al. 
2008). Assuming the detected line in DEIMOS is [Oil] at 
z = 1.53, HST G141 grism data imply a 2cr upper limit on 
[OIII]/[OII] ^0.3, which is highly unlikely for a star forming 
galaxy (Schmidt et al., 2015). Therefore, we conclude that 
the HST grism data also strongly support the z = 6.76 Lya 
interpretation of this emission line. 

We perform line fiux measurements from the DEIMOS data 
obtained during the 2013 run (with the 8301/mm grating) fol¬ 
lowing the procedure outlined in Section 2.4 of Lemaux et al. 
(2009). In short, we place two filters of width 20A on both 
sides of the emission line that are free of spectral features and 
sky line residuals to measure the background, and a central fil¬ 
ter encompassing the emission line to measure the integrated 
fiux. The background, which is fit to a linear function, is then 
subtracted from the integrated line fiux. We perform spec- 
trophotometric calibration of the DEIMOS data using two 
other compact sources (with half-light radii ~ 0."3 as mea- 
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Fig. 7.— Reduced two-dimensional (top and middle panels) and one¬ 
dimensional (bottom panel) spectrum of MACS0454-1251. The top panel 
shows the data taken on Nov 27, 2014, and the middle panel shows the data 
taken on Nov 28, 2014; we observed with the 1200 1/mm grating on both 
nights. The one-dimensional spectrum (bottom panel) was extracted from 
the data in the first night (top panel). We also plot the RMS spectra in dashed 
lines and mark the emission lines redshift if the line is to be Lya. The flux 
density values given on the ordinate are calculated in the rest-frame assuming 
the line to be Lya. 

sured from the F775W image) on the same mask with contin¬ 
uum detection in the following manner: for each object, the 
combination of slit loss and loss due to clouds was determined 
by calculating a spectral magnitude, done by correcting each 
DEIMOS spectrum for spectral response and atmospheric ex¬ 
tinction and convolving the F775W filter curve with the re¬ 
sulting spectra, and comparing this value with the magnitude 
measured in the HST image. The ratio of the fiux densities for 
each of the two sources is calculated and averaged to estimate 
the total spectral loss for this mask, which is then applied to 
RXJ1347-1216 assuming a similar half-light radius for this 
object. The reason for this assumption is, though the half- 
light radius of RXJ1347-1216 is smaller (0."2), which would, 
in principle, mean less slit loss, the size of the Lya nebula 
is known to far exceed the size of the UV continuum region 
(e.g., Wisotzki et al. 2015). For the two sources on the mask 
with which we performed the fiux calibration, the total mea¬ 
sured throughput of the slit was ^40%, lower than the ^60% 


expected for sources of this size (if they are symmetric), sug¬ 
gesting at least some departure from a photometric night. 

Using the above procedure, we measure the line fiux from 
the 8301/mm data to be 7.8 ± 0.7 x 10“^^ erg cm“^ s“^, which 
translates into a Lya luminosity of 4.1 ±0.4 x 10^^ erg s“^. 
We do not detect continuum in the spectrum, so we estimate 
the rest-frame Lya equivalent width using the object’s broad¬ 
band magnitude in F105W (on the red side of Lya) to be 
26 ± 4 A. The equivalent width uncertainty include the Pois¬ 
son noise in the central filter encompassing the sky line resid¬ 
ual, the uncertainty in the continuum, and the uncertainty in 
DEIMOS absolute flux calibration. 

We note that our measured Lya line flux from DEIMOS 
data is roughly a factor of 3 lower than that measured from 
the HST grism data, which is 2.6 ±0.5 x 10“^^ erg cm“^ s“^ 
(Schmidt et al., 2015, ApJ, in press). The difference in the line 
flux measurements could be due to the following factors: ( 1 ) 
our DEIMOS data are shallower than the HST grism data (^ 2 
hours of total integration time in the 2013 mask, versus ^ 5 
orbits of HST integration time in G102 grism), with the dif¬ 
ference in depths leading to differences in the spatial and ex¬ 
tent of the emission detected above the noise in each observa¬ 
tion, which impacts the total integrated line flux; ( 2 ) we might 
still underestimate the Lya slit loss because the true extent of 
the Lya-emitting region is much larger than the continuum- 
emitting region, while the HST slitless grism recovers more 
of the Lyo; flux; (3) there is a sky line coincident with the peak 
wavelength of the Lya emission in the DEIMOS data which 
appears slightly over-subtracted that serves to slightly lower 
the measured line flux; and (4) there could also be issues with 
contamination subtraction in the HST grism data, although 
it is unclear which direction it biases line flux measurement. 
We do expect the ground-based line flux measurement to be a 
lower limit to the space-based measurement, which is consis¬ 
tent with the measured values from DEIMOS and from HST 
grism data. 

We also measure the inverse line asymmetry \/ax as de¬ 
fined in Lemaux et al. (2009), and the line’s inverse asymme¬ 
try value {\ / ax = 0.22) is well within the range {1/ax <0.75) 
typical for convincing Lya emission. 

Using Lyo; line flux, one can also infer a star formation rate 
following Lemaux et al. (2009). The inferred star formation 
rate is strictly a lower limit, because our conversion assumes 
no attenuation of Tya photons by dust or neutral hydrogen. 
The inferred star formation rate is 1.6 ± 0.1 Mq yr“^ consis¬ 
tent with being a lower limit to the value we derive from SED 
fitting in Section 5, 17.0 ± 0.5M 0 yr“^. The Lya-inferred star 
formation rate roughly yields a Lya escape fraction of ^ 10 %. 

4.2. MACS0454-1251 

We observed MACS0454 with DEIMOS on November 28 
and 29, 2014 using the 12001/mm grating on both nights. The 
total exposure time for this mask is 7200s. We also reduce the 
DEIMOS data and extract one-dimensional spectrum follow¬ 
ing the procedure in Lemaux et al. (2009), in the same way as 
RXJ1347-1216. 

We show the reduced two-dimensional spectra of 
MACS0454-1251 in Eigure 7, and an extended emis¬ 
sion feature is clearly detected on both nights. However, 
a bright sky line residual cuts through the middle of the 
emission feature in the spectral direction and makes the line 
interpretation ambiguous. Given the width of the emission 
line, it could be either Lya at z = 6.32 or [Oil] at z = 1.39, 
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but the sky line residual makes it difficult to either confirm 
or rule out Lya or [Oil] based on line shape alone. However, 
as we will show in Section 5.2, this line is more likely to 
be Lya than [Oil] because its HST fluxes (and upper limits) 
are much better fit by a galaxy template at z = 6.32 than at 
z = 1.39 (Section 5.2), and its photometric redshift probability 
density function P{z) has very a low probability at z = 1.39 
(see the P{z) curve in Figure 10). We measure the line fluxes 
from both nights to be 1.2 ±0.2 x 10“^^ erg cm“^ s“^ (first 
night) and 8.0 ± 1.5 x 10“^^ erg cm“^ s“^ (second night), and 
these translate to Lya luminosities of 5.5 ± 0.9 x 10"^^ erg s ^ 
(first night) and 3.6 ±0.3 x 10"^^ erg s“^ (second night). From 
the line fluxes, we estimate its rest-frame equivalent widths 
(assuming Lya) using the continuum on the red side of Lya 
(estimated from its FllOW flux density) from both nights to 
be 8.2 ± 1.4 A and 5.4 ± 1.0 A. We also infer star formation 
rate lower limits to be 2.3 ± 0.4 Mq yr“^ (first night) and 
1.5 ±0.3 Mq yr“^ (second night) based on Lya fluxes, 
also fully consistent with being lower limits to the value 
derived from SED fitting in Section 5, 17.O;^4^j^M0yr“^ The 
Lya-inferred star formation rate also roughly corresponds to 
an escape fraction of ^ 10% for Lya photons. 

The inverse asymmetry value measured for the one¬ 
dimensional emission feature is ^0.5 for both masks, consis¬ 
tent with the values typical for Lya. However, the line asym¬ 
metry estimate is less reliable than that of RXJ1347-1216 be¬ 
cause we have to mask the over-subtracted skyline near the 
central wavelength of the emission feature. Based on the ob¬ 
ject’s photometric information and line asymmetry measure¬ 
ments, we identify this source as a Lya emitter at z = 6.32, al¬ 
though we are less confident with this Lya interpretation than 
we are with RXJ1347-1216. MACS0454 is not in the GLASS 
sample, so we are unable to cross check our measurements 
with HST grism data. 

5. REDSHIFT AND SPECTRAL ENERGY DISTRIBUTION 
MODELING 

In this section, we present our stellar population modeling 
of the IRAC-detected, z ^ 6 galaxy candidates using broad¬ 
band photometry. We present the modeling procedure in Sec¬ 
tion 5.1, the modeling results in Section 5.2, and we discuss 
the sources of bias and uncertainty in Section 5.3. 

5.1. Modeling Procedure 

For photometric redshift and stellar population modeling, 
we use the photometric redshift code EAZY (Brammer et al. 
2008) with stellar population templates from Bruzual & Char¬ 
iot (2003, BC03) with Chabrier (2003) initial mass function 
(IMF) between 0.1 and 100 Mq and a metallicity of 0.2 Zq. 
There is very little direct observational evidence for galaxy 
metallicity at z > 3, but limited results so far suggest that the 
majority of them have sub-solar metallicity (Maiolino et al. 
2008); we choose 0.2 Z0 for easy comparison with other 
works. The galaxy templates are generated assuming an expo¬ 
nentially declining star formation history with ^-folding time 
r ranging from 0.1 and 30 Gyr, and ages of the stellar popu¬ 
lation range from 10 Myr to 13 Gyr. For each combination of 
age and r, we implement galaxy internal dust attenuation us¬ 
ing the Calzetti et al. (2000) prescription, with the reddening 
parameter E{B-V) ranging from 0 to 1 mag to include poten¬ 
tial low-z dusty galaxy solutions. The E{B-V) grid we use 
have a step size of /SE{B-V) = 0.02 mag from E{B-V) = 0 
to 0.5 mag and a step size of /SE{B-V) = 0.1 mag from 


E{B-y) = 0.5 to 1 mag. We also use the stellar templates 
from the photometric code Le Phare (Ilbert et al. 2006)^^ in 
the fitting to check if stellar templates provide significantly 
better fits. 

Recent studies have shown that for some galaxy candi¬ 
dates, strong nebular emission lines contribute significantly to 
broadband fiuxes and therefore infiuence the inferred galaxy 
properties (e.g., Schaerer & De Barros 2010; Smit et al. 2014). 
Therefore, we use galaxy templates that include nebular emis¬ 
sion lines in the modeling. In order to calculate the expected 
line fiuxes for a given BC03 galaxy template, we calculate 
the integrated Lyman continuum flux (before dust attenua¬ 
tion) and use the relation from Leitherer & Heckman (1995) 
to calculate the expected fiuxes from hydrogen recombina¬ 
tion lines (mainly Ha, Pa/3, and Br 7 ) while assuming 
the Lyman continuum escape fraction to be zero. Non-zero 
Lyman continuum escape fraction will reduce the strength of 
optical nebular emission lines (see Inoue 2011 and Salmon 
et al. 2015). We then use the tabulated line ratios between 
H/3 and the metal lines from Anders & Alvensleben (2003) to 
calculate the metal line fiuxes for a metallicity of 0.2 Zq. For 
templates with dust attenuation, we include the dust attenua¬ 
tion effects after adding nebular emission lines. The resultant 
equivalent widths as a function of galaxy age, for r = 100 
Myr and EiB-V) = {), are shown in Figure 8, in agreement 
with Leitherer & Heckman (1995). 

In addition to nebular emission lines, we also include nebu¬ 
lar continuum emission that account for the bound-free emis¬ 
sion of HI and Hel as well as the two-photon emission of 
hydrogen from the 2s level. We follow the prescription in 
&ueger et al. (1995) (their equations 7 and 8) to calculate 
the nebular continuum flux as a function of Lyman continuum 
photon density, and we calculate the emission coefficients 
using the methods in Brown & Mathews (1970) and Nuss- 
baumer & Schmutz (1984)^^. Nebular continuum emission 
could be an important component for very young (^ lOMyr) 
starbursts and can contribute up to ^ 1/3 of the total con¬ 
tinuum just blueward of the rest-frame 4000A break. Nebu¬ 
lar continuum emission also makes the rest-frame UV slope 
redder than expected from stars alone (Schaerer & De Barros 
2010 ). 

We do not include Lya in our galaxy templates. Strong 
Lya emission could affect the LBG color selection by chang¬ 
ing the rest-frame UV broadband colors and could affect the 
derived physical properties from SED fitting (Schaerer et al. 
2011; De Barros et al. 2014). But Lya photons suffer from 
complicated radiative transfer processes and does not show 
tight correlations with the stellar population properties, so for 
simplicity we do not include Lya emission in our modeling. 

Strong gravitational lensing boosts galaxy fiuxes and in¬ 
creases the apparent SFR and stellar mass. To calculate the 
unlensed SFR and stellar mass, we use the magnification fac¬ 
tor /ibest estimated from the cluster mass models for each 
galaxy candidate at its redshift. We generate our own mod¬ 
els for MACS1149, MACS0717, MACS0454, RXJ1347, and 
RCS2327, following the procedures outlined in Bradac et al. 
(2005, 2009). In short, we constrain the gravitational poten¬ 
tial on a mesh grid within a galaxy cluster field via yf min- 

The fit using stellar templates is still done using EAZY, only the tem¬ 
plates are from Le Phare. 

The fitting formula in Nussbaumer & Schmutz (1984) is crucial to cal¬ 
culate the two-photon continuum emission between rest-frame 1216A and 
2431A, where two-photon emission dominates the nebular continuum. 
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Fig. 8.— Equivalent widths of Ua , H/3, [Oil] AA3726, 3729, and [OIII] 
A5007 that we add to the BC03 models as a function of stellar population 
age, assuming that all Lyman continuum photons are converted into nebular 
emission. The BC03 galaxy templates used in this figure have a metallicity 
of 0.2Zq , no dust attenuation, and a star-formation rate ^-folding time of 100 
Myr. Some galaxy candidates (e.g., RXJl 347-1216) require strong nebular 
emission lines to explain their observed IRAC colors, 
imization, and we adaptively use denser pixel grids near the 
core(s) of the cluster and around multiple images. We find 
the minimum values by iteratively solving a set of lin¬ 
earized equations that satisfy dx^/d^pk = 0. where is the 
gravitational potential in the kth dimension. We then pro¬ 
duce the magnification (/i) map from the best-fit gravitational 
potential map. For the rest of the clusters in our sample 
(MACS1423, MACS2129, and MACS0744), we use the pub¬ 
lic PIEMD-eNFW^^ models by Zitrin et al. (201 The 
magnification factors (and their errors) are estimated at the 
galaxy candidate positions from the z = 9 magnification maps 
except for MACS 1423-587, MACS 1423-774, MACS 1423- 
2248, and RXJ1347-1800 (which have Zbest ^ 1 so we esti¬ 
mate their //best from the z = 1 magnification map.) 

5.2. Modeling Results 

We list the best-fit galaxy properties in Table 4 and show 
the best-fit templates and the photometric redshift probabil¬ 
ity density function P{z) (while allowing redshift to fioat) in 
Figures 9 and 10. For each galaxy candidate, we estimate the 
statistical uncertainties of stellar population properties using 
Monte Carlo simulations: we perturb the photometry within 
the errors (assuming Gaussian flux errors), re-fit with the same 
set of galaxy templates, and collected the distributions of each 
best-fit property. We only perturb the fluxes where S/N > 1. 
For upper limits, we do not perturb the fluxes in our simula¬ 
tions. The systematic errors related to assumptions in initial 
mass function, galaxy metallicity, and the functional form of 
star formation history are not represented by the error bars. 
We show the distributions from Monte Carlo simulations for 
stellar mass, star formation rate (SFR), and stellar population 
age in Appendix A. From these distributions, we derive the 
confidence intervals that bracket 68% of the total probability 
in Table 4. The error bars do not include uncertainties in /ibest- 

We also show the best-fit galaxy properties for RXJl347- 
1216 and MACS0454-1251, the two galaxies that we have 

PIEMD-eNFW models use pseudo-isothermal elliptical mass distribu¬ 
tions for galaxies and elliptical NEW profiles for dark matter. 

These models are made public as high-end science products of the 
CLASH program; http : //archive . stsci . edu/prepds/clash/. 


line detections from DEIMOS data (see Section 4), when we 
fix their redshifts at their Lyct redshifts in Figure 11 (assuming 
both lines are Lyct). For RXJl347-1216, its photometric red¬ 
shift is already sharply peaked at z = 6.7, so the best-fit tem¬ 
plate and physical properties do not change after fixing its red¬ 
shift at the Lya redshift. On the other hand, MACS0454-1251 
has slightly different best-fit photometric redshift and Lya 
redshift, and we also list its best-fit properties at zcya = 6.32 
in Table 4. We also show the best-fit template at z = 1.39 in 
Figure 11 if the detected emission line is [Oil] instead of Lya 
and see that the z = 1.39 solution has a higher xj, (xj = 4.02) 
than the z = 6.32 solution (x^ = 2.90). The likelihood ratio of 

these two fits, calculated using the total x^ as e~^^ , suggests 
that the z = 6.32 model is ^ 8000 times more likely than the 
z = 1.36 model. 

The best-fit stellar masses in our sample range from 2 x 
10^ Mq to 2.9 X 10^ Mq after correcting for magnification by 
the foreground clusters and excluding the z ^ 1 interlopers. 
The stellar masses inferred from SED fitting have smaller 
statistical errors when HST fluxes are combined with IRAC 
fluxes because of the constraints on rest-frame optical emis¬ 
sion from IRAC. We show the range of stellar mass from our 
Monte Carlo simulations in the Appendix (Figure 15), and we 
find that including IRAC fluxes tighten up the possible range 
of stellar mass for every object. We also see that the range 
of stellar mass spanned by our IRAC-detected sample are not 
necessarily at the high-mass end of the observed (Gonzalez 
et al. 2011) and simulated (Katsianis et al. 2015) stellar mass 
functions at z ^ 6. In fact, the IRAC [3.6] - [4.5] color for 
several of our galaxy candidates (e.g., RXJl347-1216) sug¬ 
gest extremely young stellar population ages (^ lOMyr) and 
large equivalent widths from nebular emission lines. For these 
sources, stellar continuum emission might not dominate the 
observed IRAC fluxes, hence their true stellar masses depend 
sensitively on the equivalent widths of nebular emission lines. 
This demonstrates the combined power of strong gravitational 
lensing and deep IRAC images that allows one to measure the 
stellar mass of galaxies further down the stellar mass 
function. 

On the other hand, the SFRs and stellar population ages are 
not necessarily well constrained by SED fitting even when 
IRAC fluxes are included (see Figures 16 and 17 in the Ap¬ 
pendix). The SFRs of high-z galaxies are often calculated 
from their rest-frame UV fluxes (after correcting for dust at¬ 
tenuation), and these are often the only constraints available 
from observations. However, the UV-derived SFR depends 
critically on the amount of dust attenuation inside each galaxy, 
and the effect of dust on the rest-frame UV color is degenerate 
with the effect of stellar population age. Furthermore, UV- 
derived SFRs probe the star formation activity over the past 
^100 Myr, so it could underestimate the instantaneous SFR 
if the stellar population is younger than ^ lOOMyr; for these 
systems, nebular emission line fluxes (e.g., Ha or [Oil]) are 
better proxies for SFRs (Kennicutt 1998). We consider SFRs 
and stellar population ages more poorly constrained compared 
with stellar mass, and we will discuss the degeneracies in SED 
fitting in Section 5.3. 

IRAC fluxes also reveal four z ^ 1 interlopers from our sam¬ 
ple — MACS 1423-587, MACS 1423-774, MACS 1423-2248, 
and RXJ1347-1800 — as shown in the three bottom-right pan¬ 
els in Figure 10. All four sources have significant integrated 
probabilities at z > 6 when only HST photometry is used in 
the fitting, but the addition of their IRAC fluxes pushes their 
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Fig. 9.— Best-fit SEDs for each galaxy candidate with IRAC detections. The best-fit SEDs using only HST photometry are shown in dashed blue lines, while 
the SEDs using combined HST and IRAC photometry are shown in solid red lines. The best-fit stellar templates (fixed at z = 0) are shown in thin dotted lines. 
The photometric redshift probability density functions P(z) are shown as insets. The photometric redshifts decreases from top to bottom first, then from the left 
column to the right column. 
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TABLE 4 

Photometric Redshift and Stellar Population Modeling Results 


Object ID 

^best^ 

/4best^ 

-^stellar ^ fjj. 

(10® Mo) 

^ SFRx/^^ 
(M© yr-i) 

Age^i 

(Myr) 

sSFR® 

(Gyr-‘) 

(mag) 

Mi6oo-2.51og(/4,)k 

(mag) 




F125W-dropouts (z ~ 

9; M*y = -20.63 ± 0.36 at z - 

' 8 from Bouwens et al. 2015) 


MACS 1149-JD 

Q Q+0.1 

5 5+9-3 
‘^•‘^-0.3 

0 9+0.5 
^•^-0.4 

1 q+0.6 

1 -^-0.5 

200!go 

2 1+2-1 
^•1-0.7 

0.00 

-19.9±0.1 



F105W-dropouts (z ~ 

8; M*y = -20.63 ± 0.36 at z - 

' 8 from Bouwens et al. 2015) 


RXJ1347-1080 

MACS1423-1384 

7 3+9-3 

11m 

'^•^-0.1 

e.o+Oj 

4-Otli 

l-l-aa 

0 5+1-9 

^-^-0.5 

115 6+16-1 

290+350 
^^^-260 
< 130 

0-9^1 

lOS.O^Ogi^ 

0.00 

0.38 -0.5!^;^ 

-18.9±0.2 
-19.4 ±0.2 



F850LP-dropouts {z^l\ = - 

-20.87 ± 0.26 from Bouwens et al. 2015) 


MACS 1423-1494 

MACS0744-2088 
RXJ1347-1216g 
MACS 1423-2097 
Bullet-3^ 

MACS 1423-587 

RXJ1347-1800 

MACS 1423-774 

MACS 1423-2248 

7 1+0.3 

6.76 

6.8^;‘ 

6.8^:1 

0-l± 

0.8™| 

1 2+5A 

-1.0 

1 7+9-1 

5^8:1 

50^:1 

35^8:^ 

■^•'^-0.1 

1 9+4.0 

15^3 

-Q-1 

4 1+0.1 

12^8:1 

2^:1 

-0.1 

0 2+91 
2978:1 

^•^-2.5 

2.0+9-6 

0 1+9-2 
-Q-1 

0 1+9-1 

03^8:1 

^•^-1.0 

22 1+4-2 
^^•2-19.3 
34 0+4-9 

17 0+2-° 

5 ; 28;8 

-0.1 

1 3+1-4 

1-2-0.6 

<0.1 

< 1.2 

<20.6 

<5.2 

< 10 
20^189 
< 10 

630!‘« 

11500+12^9 
2600+||0 

1430+180 
5000^« 

J0’0'0’_4980 

45.74P 

lOS.O^j 

0 7+9-2 
'-0.6 
<0.1 

<36.0 

< 100.1 
<52.9 

0.14 -2.1!}-4 

0.16 -0.8!9:7 

0.20 -2.5!?-7 

0.00 -0.6!9:6 

0.00 

0.70 

0.00 

0.06 

0.00 

-20.1 ±0.2 
-20.7 ±0.1 
-19.3±0.1 
-19.7±0.1 
-18.9 ±0.4 


RCS2327-1282 

MACS0454-1817 

MACS0454-1251 

MACS0454-125lj 


F814W-dropouts (z ~ 6-7; = -20.87 ± 0.26 from Bouwens et al. 2015) 


1 q+0.4 
- 0.1 

2.1 


0.5; 






19.0; 


-5.8 

,+7.9 

-5.4 


400!}oo 

< 30 

30!g 


o 9+1.6 

105 . 0 + 0 ,' 

O 

■^^•'^- 35.3 


0.00 

0.28 

0.12 

0.08 


- 3 . 0 ! 0 - 

-1.4; 

- 1.6 

-1.6; 


;8:i 


-20.9 ±0.1 
-19.3±0.1 
- 20.8 ± 0.1 
-20.9 ±0.1 




F775W-dropouts (z ^ 

6; M^y = -20.94 ± 0.20 from Bouwens et al. 2015) 



MACS 1149-274 

MACS 1149-1204 

j'li 

2*-±0.1 

1 o+O.l 1 O+0.3 

i.O-Oi ^-^-0.8 

7.3!3'^ 0.08 

12.9!|7 sO^^o 10.2!!” 0.08 

-1 6+9-1 
_1 578:1 

^•2’-0.5 

-21.2±0.1 
-20.7 ±0.2 


^ Zbest is the photometric redshift using BC03 galaxy templates except for RXJ1347-1216 and MACS0454-1251, for which we identify Lya emission at Zspec = 6.76 and 
Zspec = 6.32, respectively. 

^ Lensing magnification factor estimated from the galaxy cluster mass models mentioned in Section 5.1. For MACS 1423-587, MACS 1423-774, MACS 1423-2248, and 
RXJ1347-1800, we estimate their ^tbest from the magnification map at z = 1. 

The intrinsic stellar mass and SFR assuming fi = fihest- To use a different magnification factor /i,, simply use = [i/ /ibest, where /ibest is the best magnification factor 
we adopt for each object. When the best-fit SFR is zero, we report the 68% upper limit. 

^ Time since the onset of star formation. For the sources with best-fit age equal to 10 Myr, the youngest template allowed in our models, we report the 68% upper limit of 
the age from Monte Carlo simulations. 

^ Specific star formation rate = SFR / stellar mass. When the best-fit sSFR is zero, we report the 68% upper limit. 

® All fits are performed at z = 6.76, the Lya redshift. The confidence intervals reflect the maximal range returned from the simulations because the distributions for this 
object are highly skewed. 

^ The best-fit color excess E(B-V) of the stellar emission from our SED modeling. The dust attenuation at rest-frame 1600A can be calculated using the dust attenuation 
curve from Calzetti et al. (2000) as Aiboo = 9.97 x E{B-V). 

^ Measured rest-frame UV slope /3 fromZ/^TAYFCS broadband fluxes, assuming the candidates are at z = Zbest- Wedo not measure /3 for MACS 1423-587, MACS1423-774, 
MACS 1423-2248, and RXJ1347-1800 because they have Zbest ~ L We also do not measure /3 for MACS1149-JD because it has only 1 filter (F160W) that samples the 
rest-frame UV continuum. 

J All fits are performed at z = 6.32, its Lya redshift. 

^ Rest-frame 1600 A absolute magnitude assuming z = Zbest and fi = fihsst- 
^ First reported by Ryan et al. (2014); included here for completeness. 
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Fig. 10.— Same as Figure 9, for the remaining eight IRAC-detected z > 6 galaxy candidates. The two galaxy candidates with Zbest ~ 1 are shown here. 

























































































































16 



the rest-frame 4000A break instead of the hya break. This 
demonstrates the value of IRAC detections in discriminating 
between genuine z ^ 6 galaxies and lower-z interlopers. 

In Figure 9 and 10 we also show the best fit stellar templates 
from Le Phare. Most of the sources are better fit by galaxy 
templates than by stellar templates, although there is only one 
case, RXJ1347-1800, where both templates provide similarly 
good fits (xj, = 0.57 for galaxy templates and x?. = 0.62 for 
stellar templates). For all the galaxy candidates, their best- 
fit stellar templates are either brown dwarfs or low-mass stars 
from Chabrier & Baraffe (2000). We also check if our sample 
contains X-ray detected sources that could have significant 
contributions from AGNs and we do not find evidence that 
any of our galaxy candidates have AGNs. However, low-level 
AGN activities are still possible, and the stellar population 
properties we infer from SED modeling will depend on how 
much (if any) AGN contribution there is to their broadband 
fiuxes. 

5.3. Modeling Biases and Uncertainties 

The biases and uncertainties of SED fitting are well doc¬ 
umented in the literature (e.g., Papovich et al. 2001; Lee 
et al. 2009), and most sources of systematic errors come from 
model assumptions. In general, stellar mass has the small¬ 
est systematic errors (~ 0.3 dex), but uncertainties in galaxy 
star formation history can lead to large biases in galaxy age 
and star formation rate (Lee et al. 2009). In additional to the 
usual culprits of systematic errors (e.g., star formation his¬ 
tory, initial mass function), another important source of sys¬ 
tematic error is the nebular emission line ratios. We find that 
nebular emission line contributions to broadband fiuxes are 
important for a subset of our sample, but direct observational 
constraints of rest-frame optical nebular emission line ratios 
of z > 6 galaxies will not arrive until the launch of JWST. Un¬ 
certainties in the amount of dust attenuation also complicates 
the interpretation of the best-fit parameters, as we demonstrate 
below. 

We use MACS1149-JD to show how uncertainties in dust 
attenuation can lead to uncertainties in star formation rate in 
Figure 12. We estimate the number density contours from 
our Monte Carlo simulations (with 1000 realizations) when 


Fig. 12. — Distribution of galaxy age v.s. star formation rate for 
MACS 1149-JD (Zphot = 9.3) from our Monte Carlo simulations. We plot the 
estimated number density contours for Monte Carlo realizations with differ¬ 
ent ranges of E(B-V) values in the central panel to show the correlation 
between galaxy age, dust attenuation, and star formation rate. We also show 
the marginalized histograms of star formation rate on the right, and we show 
the marginalized histograms of galaxy age on the top. The star formation rate 
is the intrinsic value assuming a magnification factor of 5.5 for MACS1149- 
JD. 

IRAC fiuxes are included in the modeling. We show the 
number density contour of each E{B-V) value in the cen¬ 
tral panel and the marginalized histograms for galaxy age (on 
top) and star formation rate (on the right). The star forma¬ 
tion rate histogram shows a single peak at \.9Mq yr“^ assum¬ 
ing a magnification factor of 5.5, but there is a long tail to 
higher star formation rates that extends one order of magni¬ 
tude. The long tail in star formation rate corresponds to higher 
dust attenuation templates {E{B-V)> \ green dashed con¬ 

tours in the central panel) as opposed to the peak of the his¬ 
togram {E{B-V) < 0.1; solid blue contours in the central 
panel). We note that because of the high photometric red- 
shift of MACS1149-JD (Zphot = 9.3), its UV continuum be¬ 
tween rest-frame 1250 and 2600 A is not well sampled by 
HST and IRAC photometry; the addition of X-band photome¬ 
try should help constrain the amount of dust attenuation inside 
this galaxy. 

As a model-independent check on the inferred dust attenu¬ 
ation, we measure the rest-frame UV slope /3 for each galaxy 
candidate and list the values in Table 4. We can only measure 
(3 when a galaxy candidate has at least two filters sampling the 
UV continuum (between rest-frame 1250 and 2600A); there¬ 
fore, we do not measure p for MACS1149-JD (which only 
has F160W that samples UV continuum) and for MACS1423- 
587, RXJ1347-1800, MACS1423-774, and MACS1423-2248 
(which have photometric redshifts ^1). We measure (3 by us¬ 
ing a power-law spectrum fx oc convolving the spectrum 
with the filter curves that sample the UV continuum, and find¬ 
ing the (3 that best matches the observed fiuxes. The uncer¬ 
tainties are quantified in bootstrap Monte Carlo simulations, 
and we show the E{B-V) values from SED fitting v.s. p in 
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Fig. 13. — Best-fit E(B-V) v.s. (3 for the IR AC-detected 6 < z < 10 sam¬ 
ple. We also show the expected /3 from a given E(B-V) value assuming a 
dust attenuation law from Calzetti et al. (2000) and the empirical relations 
between A 1600 to (3 from Meurer et al. (1999) (for solar metallicity) and from 
Castellano et al. (2014) (for sub-solar metallicity). We randornly shift the 
best-fit E(B-V) for each galaxy candidate around its best-fit value by no 
more than 0.01 for clarify. As a whole sample, the measured (3 values are 
consistent with the expected (3 values from E(B-V), although the scatter is 
large and the dust attenuation for each galaxy candidate is poorly constrained. 


Figure 13. 

In Figure 13, we also show the expected values of P 
given values of E{B-V) using two different empirical cal¬ 
ibrations. To calculate the expected P, we use the Calzetti 
et al. (2000) dust attenuation law to calculate the amount of 
dust attenuation at rest-frame 1600A from E{B-V)\ Ai6oo = 
9.97 xE{B-V). Then we use the relation between Ai6oo and 
/3 from Meurer et al. (1999) (for solar metallicity; Ai6oo = 
4.43 + 1 .99/3) and Castellano et al. (2014) (for sub-solar metal¬ 
licity; Ai6oo = 5.32;^q 37 + 1.99y5) to calculate the expected p. 
We show the Meurer et al. (1999) relation as a solid line and 
the Castellano et al. (2014) relation as a dashed line in Figure 
13. The scatter of the measured /3 of our sample is larger than 
the difference between the two empirical calibrations, but the 
distribution is roughly consistent with both calibrations. The 
agreement means that the E(B-V) values derived from SED 
fitting is not strongly biased on average, but the large scatter 
also suggests that for each individual galaxy, the dust attenu¬ 
ation is still poorly constrained. 

6. IRAC COLORS AND STRONG NEBULAR EMISSION LINES 

Recent works suggest that at least in a subset of high-z 
galaxies, strong nebular emission lines (most notably Ha, 
up, [OIII] A5007, and [Oil] A3727) with rest-frame equiva¬ 
lent widths ^ 200 A or higher contribute significantly to their 
IRAC fluxes (e.g.. Shim et al. 2011; Schaerer & De Barros 
2009; De Barros et al. 2014; Smit et al. 2014). The galaxies 
with extreme nebular emission line strengths are most likely 
starbursts younger than 100 Myr; such galaxies are also being 
found in increasing numbers at z ^ 2-3 (e.g., van der Wei et al. 
2011; Atek et al. 2011). If a large number of such galaxies 
exist at z ^ 6, strong nebular emission lines in the rest-frame 
UV/optical wavelengths need to be included in the stellar pop¬ 
ulation modeling. 

Within certain redshift ranges, unusual IRAC [3.6]-[4.5] 
colors can be tell-tale signs of strong nebular emission lines. 
Shim et al. (2011) identified 47 galaxies at z ^ 4 that have 
bluer [3.6] - [4.5] colors than those expected from stellar con- 
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Fig. 14.— IRAC [3.6]-[4.5] color as a function of redshift calculated 
from BC03 model spectra. We show the [3.6] - [4.5] colors for 0.2 Z0 Stellar 
population models that are 500 Myr old without nebular emission lines (thin 
dot-dashed curve), 10 Myr old without emission lines (thin dotted curve), 
and a 10 Myr old model with nebular emission lines (thick solid curve). In 
our implementation, the dust-free 10 Myr model with nebular emission lines 
has Ha, ¥1(3, and [OIII]A5007 equivalent widths of 1087A, 182A, and 868A, 
respectively; for the 10 Myr model, the equivalent widths change by < 3% 
within the range of ^-folding time r that we adopt. In addition to the fidu¬ 
cial 0.2 Zq models, we also show the expected colors of a 0.02 Zq, 10 Myr 
old model with nebular emission lines for comparison (thick dashed curve); 
the more metal-poor model seems to reproduce the colors of RXJ1347-1216. 
The measured [3.6] - [4.5] colors of two sources in our sample that likely 
have strong emission lines are shown in stars, and the rest of the sample is 
shown in black circles. We show the known Lya emitters (LAEs) infilled 
symbols, including the two LAEs reported in this work (see Section 4). We 
also show other published z > 6 LAEs with measured [3.6] - [4.5] colors — 
z8_GND_5296 (Finkelstein et al. 2013, F13), HCM6A (Chary et al. 2005, 
C05), GN-108036 (Ono et al. 2012, 012), CR7 (Sobral et al. 2015, S15), 
EGS-zs8-l (Oesch et al. 2015, 015), EGS-zs8-2 (Roberts-Borsani et al. 2015, 
RB15), and EGSY-2008532660 (Zitrin et al. 2015, Z15) — in filled sym¬ 
bols. The [3.6]-[4.5] color of z8_GND_5296, HCM6A, EGS-zs-8-2, and 
MACS 1423-1494 are hard to reproduce by the stellar population models that 
we adopt, but they come close to the expected colors of a Type 2 (obscured) 
AGN template from (Polletta et al. 2006, P06). 


tinuum alone, and they categorized these galaxies as Ha emit¬ 
ters because the SED bumps at 3.6/im are likely due to strong 
Ha emission. More recently, Smit et al. (2014) and Smit et al. 
(2015) presented z ^ 6.6-7.0 galaxy candidates with unusu¬ 
ally blue [3.6]-[4.5] colors as evidence for strong contribu¬ 
tions from [OIII] and to the 3.6/im fluxes. The colors 
of these peculiar objects usually can only be reproduced with 
model SEDs that include strong nebular emission lines. 

In Eigure 14, we compare the IRAC [3.6] - [4.5] colors of 
our sample with a range of model predictions. We use our 
fiducial SED model (BC03) to generate the redshift evolution 
of [3.6] - [4.5] color at 500 Myr old without nebular emission 
lines (thin dot-dashed curve, roughly the age of the universe at 
z = 9.5), 10 Myr old without nebular emission lines (thin dot¬ 
ted curve), and 10 Myr old with nebular emission lines (thick 
solid curve). The 10 Myr old model with nebular emission 
lines have equivalent widths 1087A, 182A, and 868A for Ha, 
H/3, and [OIII]AA4959,5007, respectively. Here we assume 
the star formation p-folding time scale r to be 100 Myr, but 
the [3.6] - [4.5] color does not change significantly when dif¬ 
ferent values of r are used. 

The most prominent feature in Eigure 14 is the “dip” in 
[3.6] - [4.5] for a 10 Myr old starburst with nebular emission 
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lines at z ^ 6.8 due to the contributions from [OIII] and H/3 

— the same feature that Smit et al. (2014) utilized to identify 
strong nebular emission line objects within 6.6 < z < 7. In our 
sample, only RXJ1347-1216 has a photometric redshift ~ 6.8 
and a very blue [3.6] - [4.5] color. This source has a best-fit 
age of 10 Myr, the youngest age included in our templates. 
Extremely young stellar populations are expected to generate 
a large number of ionizing photons, so if these sources are in¬ 
deed ^10 Myr old starbursts, they might also have high Lya 
luminosities around star forming regions. We already suc¬ 
cessfully identified one of the three sources (RXJ1347-1216) 
as a z = 6.76 Lya emitter (LAE; see Section 4); we do not 
identify other sources at z ^ 6.8 with blue [3.6] - [4.5] color 
that could also be strong line emitters in our sample. 

In Eigure 14 we also show the redshift evolution of [3.6] - 

[4.5] color for a 10 Myr old, 0.02Zq model (thick dashed 
curve), and it predicts a bluer [3.6] - [4.5] color at z ^ 6.8 
(as blue as ^ -1.4 mag) than the 10 Myr old, 0.2Zq model. 
The IRAC colors of the 0.02Zq model at z ^ 6.8 show bet¬ 
ter agreements with the three sources mentioned above than 
the O. 2 Z 0 model, which suggests that these sources might 
have lower metallicities than our fiducial model. We note that 
the nebular emission line properties of individual galaxies are 
highly uncertain (and are sensitive to metallicity), so any con¬ 
straint on metallicity is preliminary. 

Our fiducial galaxy SED models also predict that young 
starbursts with strong nebular emission lines should have red 

[3.6] - [4.5] colors at 7.0 < z < 7.5. In this redshift range, 
[Oil] and [OIII] move into IRAC chi and ch2, respectively, 
and the combined [OIII]-rH/3 line fiux is expected to be ^ 4 
times higher than [O II] in our implementation (Anders & 
Alvensleben 2003); the expected [3.6]-[4.5] color reaches 

- 0.5 mag within 7.0 < z < 7.5. MACS1423-1494 (zphot= 
7.3) has a photometric redshift and measured [3.6] - [4.5] 
color that are close to the model prediction, although its red 

[3.6] - [4.5] color is hard to reconcile even with the 10 Myr 
old galaxy model. Based on its photometric redshift and un¬ 
usually red [3.6] - [4.5] color (it is also best-fit by a 10 Myr 
old galaxy SED), we identify this source as another prime 
candidate for Lya emission. We note that Lya photons are 
subject to complicated radiation transfer effects both inside 
and outside of galaxies, so it is far from guaranteed that these 
sources will have detectable Lya emission. But they are likely 
LAE candidates (compared with other high-z galaxies) based 
on their photometric redshifts and IRAC colors. 

We also compare our galaxy model-predicted IRAC colors 
with other z ^ 6.5 LAEs with published IRAC colors in Eigure 
14. The other LAEs include HCM 6 A from Hu et al. (2002) 
(z^ = 6.56, where z^ is the spectroscopic redshift determined 
by Lya emission), CR7 from Sobral et al. (2015) (z^ = 6.60), 
GN-108036 from Ono et al. (2012) (z, = 7.21), EGS-zs8-2 
from Roberts-Borsani et al. (2015) (z^ = 7.48), z8_GND_5296 
from Einkelstein et al. (2013) (z^ = 7.51), EGS-zs 8 -l from 
Oesch et al. (2015) (z, = 7.73), and EGSY-2008532660 from 
Zitrin et al. (2015) (z^ = 8 . 68 ). All of these LAEs have IRAC 
colors that strongly suggest high nebular emission line equiv¬ 
alent widths (most likely [OIII] and at this redshift range), 
because they lie along the curve traced by a dust-free, 0.2 Zq, 
10 Myr stellar population. Eor example, Einkelstein et al. 
(2013) argued that the red IRAC color of z8_GND_5296 is 
due to the galaxy’s strong [OIII]-l-H/3 emission lines in IRAC 
ch2, and they inferred the [OIII] A5007 equivalent width to 
be 560-640 A from photometry. The IRAC colors of these 


z ^ 6.5 LAEs corroborate the recent findings that many galax¬ 
ies detected at z ^ 6 likely have high nebular emission line 
equivalent widths. 

Two notable cases among the group of LAEs in Eigure 
14 are MACS 1423-1494 and HCMbA^^. HCM6A was found 
in the vicinity of a massive galaxy cluster Abell 370^^ and 
has a measured [3.6]-[4.5] color of 1.0 ±0.4 mag, signifi¬ 
cantly redder than the [3.6] - [4.5] color predicted by a 10 Myr 
stellar population model at its redshift (z^ = 6.56). The red 

[3.6] - [4.5] color suggests a very high Ha/([OIII]-rH/3) ratio, 
which is unexpected (but not impossible) for a young, low- 
metallicity stellar population. In order to explore other pos¬ 
sibilities to explain the red [3.6] - [4.5] colors of both LAEs, 
we plot the predicted [3.6] - [4.5] colors of a Type 2 obscured 
AGN template from Polletta et al. (2007). This obscured AGN 
template includes a dust attenuation of Ay = 4 mag that fits 
the obscured AGN SW 104409 (z = 2.54; Polletta et al. 2006), 
and its color trajectory in redshift is shown as a thick dotted 
curve in Eigure 14. Interestingly, the predicted [3.6]-[4.5] 
colors of an obscured AGN agrees quite well with the colors 
of both MACS 1423-1494 and HCM6A, and z8_GND_5296 
and EGS-zs8-2 also have marginally consistent IRAC colors 
with this obscured AGN template. If these sources indeed har¬ 
bor obscured AGNs (like SW 104409), the red [3.6]-[4.5] 
colors will be primarily due to large dust attenuation in the 
rest-frame optical, while the blue rest-frame UV colors come 
from the scattered light of the central QSO emission. Ob¬ 
scured AGN is an intriguing possibility to consider for these 
sources, although so far no direct evidence exists that any of 
these sources have significant fiux contributions from an ob¬ 
scured AGN. 

To sum up, we identify three sources in our sample at 
z ^ 6.7 and z ^ 7.3 as likely young starbursts with very strong 
nebular emission lines based on their IRAC colors. We detect 
Lya emission in one of them, RXJ1347-1216, during our re¬ 
cent DEIMOS observations, and we plan to follow up all the 
other three sources for their potential Lya emission. 

7. SUMMARY 

In this work, we present the constraints on the 6 < z ^ 10, 
IRAC-detected galaxy candidates behind eight strong-lensing 
galaxy clusters from SURES UP. Six of the clusters are in the 
CLASH sample, and two are in the Hubble Erontier Eields 
sample. We summarize our findings as follows: 

• We find a total of 17 galaxy candidates using the Ly¬ 
man break color selection that have S'/A > 3 in at least 
one IRAC channel. The photometric redshifts in our 
sample range from 5.7 to 9.3, and we identify four 
galaxy candidates (MACS 1423-587, RXJ1347-1800, 
MACS 1423-774, and MACS 1423-2248) as likely z - 1 
interlopers after including their IRAC fiuxes in the SED 
modeling. We find the largest number (6) of IRAC- 
detected galaxy candidates in MACS 1423. 

• Erom our Keck spectroscopic observations, we identify 
one secure Lya emitter at z = 6.76 (RXJ1347-1216) 
and one likely Lya emitter at z = 6.32 (MACS0454- 
1251). The line equivalent widths, assuming they are 
both Lya, are 26±4A (RXJ1347-1216) and 6.8 ± 1.7A 

HCM6A was first reported by Hu et al. (2002), and later Chary et al. 
(2005) published its IRAC fluxes. 

Abell 370 is one of the Hubble Frontier Fields cluster. 
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(MACS0454-1251, averaged over two nights). We in¬ 
fer lower limits of their star formation rates from their 
hya line fluxes and find them to be consistent with the 
star formation rates from SED fitting. 

• We infer the physical properties of our sample galax¬ 
ies using Bruzual & Chariot (2003) galaxy templates 
and add nebular emission lines to the templates. Un¬ 
der our SED modeling assumptions (0.2 Zq, Chabrier 
IME, exponentially decaying star formation history, and 
nebular line emission), the stellar masses of our sam¬ 
ple range from 0.2-2.9 x 10^ Mq (excluding the three 
likely z ^ 1 interlopers) when we use the best avail¬ 
able magnification factors for each galaxy candidate. 
The magnification-corrected rest-frame 1600 A abso¬ 
lute magnitude (Mi6oo; see Table 4) of our sample 
ranges from -21.2 to -18.9 mag. The range of intrin¬ 
sic UV luminosity probed here is slightly fainter than 
the knee of UV luminosity functions at 6 < z < 10, 
which have between ^ -20.6 and ^ -21.6 mag 
(e.g., Bouwens et al. 2015; Einkelstein et al. 2015), 
showing that galaxy clusters’ strong tensing power al¬ 
lows us to start probing the more typical UV luminosi¬ 
ties. The range of intrinsic stellar mass probed here is 
also close to the knee of the stellar mass functions at 
this redshift range (e.g., Gonzalez et al. 2011; Katsia- 
nis et al. 2015). Some galaxies in our sample are best 
fit by extremely young (^10 Myr old) templates and 
others best fit by more evolved (up to ^ 700 Myr old 
at z ^ 7) templates, suggesting that the IRAC-detected 
sample contains both very young galaxies with strong 
nebular emission lines and more evolved and massive 
galaxies at 6 < z ^ 10. 

• Erom the photometric redshifts and IRAC colors, we 
identify two galaxy candidates that likely have strong 
(rest-frame optical) nebular emission lines: RXJ1347- 
1216 and MACS 1423-1494. Both sources are best 
fit by the youngest (10 Myr old) galaxy templates in¬ 
cluded in our modeling and are prime targets for spec¬ 
troscopic observations. We already identified one of 
them (RXJ1347-1216) as a Lya emitter, and we will 
target the other two in our future spectroscopic observa¬ 
tions. Other galaxies in the sample lie in the part of the 
redshift-IRAC color space that makes it hard to infer 


their nebular emission line strengths; namely, they are 
within the redshift range that both IRAC bands could 
have contribution from strong nebular emission lines 
such as [OIII], H/3, and Ha, and their IRAC colors may 
not be very different from those of pure stellar contin¬ 
uum. 

The IRAC fluxes provide important information about the 
galaxies at 6 < z ^ 10, because it is the only probe of their 
rest-frame optical emission that we have at the moment. The 
IRAC-detected galaxies may not be representative of the 
entire galaxy population at z ^ 6, but their IRAC colors do 
provide a more effective way to select spectroscopic targets 
for redshift confirmation. IRAC fiuxes and meaningful upper 
limits can also distinguish some lower-redshift galaxies from 
high-z dropouts and are important for constructing clean 
z ^ 6 galaxy samples. 
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APPENDIX 

DISTRIBUTIONS FROM MONTE CARLO SIMULATIONS 

Here we show the distributions of stellar mass, star formation rate, and stellar population age (assuming an exponentially 
declining star formation history with ^-folding time between 0.1 and 30 Gyr) in Eigures 15, 16, and 17, respectively. In all panels, 
the distributions from using HST photometry only are shown in gray filled histogram, while the distributions from combining 
HST and Spitzer photometry are shown in red histogram. 
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Fig. 15.— Stellar mass distributions derived from Monte Carlo simulations for each IRAC-detected z > 6 galaxy candidate. The stellar mass values have 
been scaled by our best estimates of magnification factors yUbest- The distributions from combining HST and IRAC photometry are shown as the unfilled red 
histograms; the distributions from HST photometry only are shown as filled gray histograms. The best-fit stellar masses from Table 4 are shown as the vertical 
dashed lines. 
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Fig. 16.— Star formation rate (SFR) distributions derived from Monte Carlo simulations for each IRAC-detected z > 6 galaxy candidate. The SFR values 
have been scaled by our best estimates of magnification factors yUbest- All SFRs below 0.01 Mq yr“^ are set to 0.01 Mq yr“^ for clarity. The distributions from 
combining HST and IRAC photometry are shown as the unfilled red histograms; the distributions from HST photometry only are shown as filled gray histograms. 
The best-fit SFR from Table 4 are shown as the vertical dashed lines. 
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Fig. 17.— Model stellar population age distributions derived from Monte Carlo simulations for each color-selected, IRAC-detected z > 6 galaxy candidate. The 
minimum age included in the template library is 10 Myr. The distributions from combining HST and IRAC photometry are shown as the unhlled red histograms; 
the distributions from HST photometry only are shown as hlled gray histograms. The best-ht age from Table 4 are shown as the vertical dashed lines. 
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